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ABSTRACT 

The thermonuclear explosion of a massive white dwarf in a Type la supernova explosion is char- 
acterized by vastly disparate spatial and temporal scales. The extreme dynamic range inherent to 
the problem prevents the use of direct numerical simulation and forces modelers to resort to subgrid 
models to describe physical processes taking place on unresolved scales. 

We consider the evolution of a model thermonuclear flame in a constant gravitational field on a 
periodic domain. The gravitational acceleration is aligned with the overall direction of the flame 
propagation, making the flame surface subject to the Ray leigh- Taylor instability. The flame evolution 
is followed through an extended initial transient phase well into the steady-state regime. The prop- 
erties of the evolution of flame surface are examined. We confirm the form of the governing equation 
of the evolution suggested by Khokhlov (1995). 

The mechanism of vorticity production and the interaction between vortices and the flame surface 
are discussed. Previously observed periodic behavior of the flame evolution is reproduced and is 
found to be caused by the turnover of the largest eddies. The characteristic timescales are found to be 
similar to the turnover time of these eddies. Relations between flame surface creation and destruction 
processes and basic characteristics of the flow are discussed. We find that the flame surface creation 
strength is associated with the Ray leigh- Taylor time scale. Also, in fully-developed turbulence, the 
flame surface destruction strength scales as l/£ 3 , where L is the turbulent driving scale. 

The results of our investigation provide support for Khokhlov's self-regulating model of turbulent 
thermonuclear flames (Khokhlov 1995). Based on these results, one can revise and extend the original 
model. The revision uses a local description of the flame surface enhancement and the evolution of the 
flame surface since the onset of turbulence, rendering it free from the assumption of an instantaneous 
steady state of turbulence. This new model can be applied to the initial transient phase of the flame 
evolution, where the self-regulation mechanism yet to be fully established. Details of this new model 
will be presented in a forthcoming paper. 
Subject headings: supernovae: general - turbulence 



1. INTRODUCTION 

Thermonuclear combustion was originally pro- 
posed as the mechan i sm behind Type la super- 
novae (|Hovle fc Fowled ll96Ct lArnettl ^969). How- 
ever, this detonation-based model was soon shown 
to be unable to explain a host of obse r vation s 
(lArnettl Il973 lOstriker Richstone. fc ThuarJ 11974ft . 
iNomoto. Sugi moto. fc Ned (|1976T ) recognized the possi- 
bility that a deflagration, a subsonic mode of combustion 
following a mild ignition of stellar material in the center 
of a massive white dwarf, might help to ameliorate these 
deficiencies. Flames characteristic of a deflagration 
are subject to the Ray leigh- Taylor instability in the 
gravitational field in the interior of a white dwarf, 
significantly complicating any detailed calculation of 
this mode of burning. 

The parameterized one-dimens ional study of 
INomoto. Thielemann, fc Yokoil (^984) indicated that a 
successful explosion model may require a substantial 
increase in the speed of the deflagrating front. Such an 
increase in burning speed was observed in the multi- 
dimensional models of Miill er fc Arnettl (|1982t [l986), 
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though Miill er fc Arnettl (|1986ft recognized that their 
models lacked a proper description of turbulence (and 
also lacked details of the microphysics) taking place on 
scales unresolved in their simulations. 

The structure and dynamics of deflagration fronts on 
Rayl eigh- Tay l or- do minated scales were studied in detail 
bv iKhokhlovl (|1995l) . In particular, Khokhlov proposed 
that when a steady state of fully-devloped turbulence is 
assumed, the evolution of the flame propagation becomes 
independent of details of the microphysics and is gov- 
erned by a self-regulating mechanism reflecting an inter- 
play between the Rayleigh- Taylor instability and nuclear 
burning. Khokhlov (1995) provides scaling relations for 
the turbulent flame speed and these scalings are used to 
construct a subgrid-scale (SGS) model which, combined 
with a flame capturing scheme, offers a powerful tool for 
studying deflagrations in Type la supernovae. 

The use of subgrid-scale models in modeling of Type la 
supernovae is necessitated by the vastly disparate tem- 
poral and spatial scales present in the problem. Without 
exception, all groups involved in studying thermonuclear 
supernovae use some sort of subgrid-scale model. In the 
Khokhlov self-regulation scenario, the turbulent flame 
speed is primarily determined by the speed of large-scale 
motions driven by the Rayleigh- Taylor instability. This 
speed is oc ^ AgL, where A is the Atwood number, g is 
the local gravitational acceleration, and L is the driving 
scale. In this subgrid model, it is assumed that a steady 
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state of fully-developed turbulence obtains instanta- 
neously, leading to a state of self regulation where the 
interplay of the RT instability and the local progress 
of the flame serves to render the turbulent flame speed 
roughly constant in time. However, the model does not 
take into account the initial transient phase nor the 
unsteadi ness of local turbulent v e locitie s. The subgrid 
model of lNiemever fc Hillebrandtl ([1995ft uses the frame- 
work of Kolmogorov turbulence theory as a starting 
point. In their model, the burning speed is locally de- 
termined from the turbulent kinetic energy, k sgs , which 
is evolved thr ough the use of a local transport equation 
(|Clementlll 993ft . This model has been used extensively 
to study thermonuclear deflagrations both without 

( Nie meve r. Hillebrandt. fc Woosle vl 1 19961) and , later, 

with (|Reinecke. Hillebrandt. fc Niemeverlll999l 12002ft a 
flame front tracking scheme. Despite the fundamental 
differences in these subgrid models, the results obtained 
(explosion energetics, timescale, and ejecta morphol- 
ogy) are quite simil ar for deflagrations in the s elf- 
regulation regime (Reinecke, Hillebrandt, & Nieme 



egime 
2002; 



Gamezo e 



t all 
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LSchmidt. Niemever. Hillebrandt. fc Ronkel 12006ft . Ac- 
cordingly, we can conclude that in the self-regulation 
regime differences between specific subgrid models do 
not lead to gross differences in the flame evolution. 

Our aim here is to improve the original Khokhlov 
subgrid-scale model (hereafter KSGS) and develop a new 
subgrid model that can be applied not only to the self- 
regulation regime but also to the initial transient phase 
before self regulation has been fully established. The 
inclusion of a description of this transient phase in the 
subgrid model is especially important for the simulation 
of the early evolution of deflagrations in Type la super- 
novae. 

To develop the new subgrid model, we have first con- 
ducted a survey of thermonuclear model flame simula- 
tions on periodic domains. These results allow us to de- 
velop a deeper understanding of the flame dynamics on 
small scales and provide a foundation for the improved 
scheme. 

2. NUMERICAL METHODS AND MODELS 

In our study of model thermonuclear flames, we con- 
sider the evolution of the reactive Eulerian equations in 
Cartesian geometry. Numerical solutions were obtained 
with the FLASH code (|Frvxell et al.ll2000ft . The majority 
of the models presented here were computed in three- 
dimensions, while a limited number of models have been 
calculated in two-dimensions. 

2.1. Numerical Methods 

We assume a constant gravitational acceleration, g Zl 
aligned with and pointing along the negative z coordi- 
nate. The computational domain is a tube (a channel in 
2D) with square cross-section. We use periodic boundary 
conditions along the lateral (i.e. (x,y)) directions. We 
impose reflecting boundary conditions at the bottom of 
the tube, while allowing outflow through its top bound- 
ary. The time steps in our simulation were limited with 
Courant factor of 0.6. 

In order to save computational resources, we use 
flash's adaptive mesh refinement (AMR) capability and 
only resolve regions actively participating in combustion. 



In these regions, we track flow structures characterized 
by large jumps in density, total velocity, or progress vari- 
able. 

At the initial time, a one-dimensional flame front struc- 
ture is interpolated onto the grid with the flame front 
position, Zf, defined as 

Zf = zo + Af cos(2ttx / X x ) cos(2ity / \ y ) , 

where zq is the unperturbed position of the front, Af 
is the amplitude of the perturbation, and = \ y = A is 
the wavelength of the perturbation in the lateral direc- 
tions. The fuel filling the tube ahead of the flame front is 
initially at rest and composed entirely of equal mass frac- 
tions of carbon and oxygen. The fuel's density is 1 x 10 8 
g/cm 3 and its temperature is 5 x lO 7 ^. We consider here 
only a one-step reaction in the burning front with a con- 
stant prescribed energy release of 7 x 10 17 erg/g (i.e. the 
total energy that would be released in taking the initial 
carbon-oxygen mix to a compo sition of all iron). We use 
a Helmholtz equation of state ([Timmes fc Swestvl feOQO) 
suitable for the degenerate conditions considered here. 

The original flash code has been equipped with a cus- 
tom implementa ti on of the flame - captu ring algorithm of 
iKhokhlovl {l995); Ga mezo et all ([2003ft . In this frame- 
work, the flame front is described as a diffuse, thick inter- 
face of a progress variable, /, whose evolution is governed 
by an advection-diffusion-reaction (ADR) equation: 



dl 
dt 



U-Vf = K\7 2 f + R, 



(1) 



Here / is constrained to lie between / = (unburned 
material, or fuel) and / = 1 (completely burned material, 
or ash). The diffusion and reaction coefficients, K and 
R, are given by 



K 
R- 



const, 



const, if f < f < 1 , 
otherwise , 



where fo = 0.3. Equation describes a reactive front 
of thickness S ~ (K/C) 1 ^ 2 propagating with the speed 
D = (KC/f ) 1/2 . With coefficients K and C defined as 



K = D(f3Ax)^ 



C 



D 



(J3Ax) 



and /3 = 1.5, the front diffuses with a prescribed con- 
stant speed D and spreads over ^3 — 4 zones. For 
the advection part of the A DR scheme we use the PPM 
(|Colella fc WoodwardlH98l solver in flash. For a more 
detailed description of the flam e-capturing method an d 
its characteristics see Khokhlo^JTHa) and lzhiglol (12005). 
Verification of the flash implementation of this scheme 
has b een presented in lyiadimirova, Weirs, <fe Hyzhik| 
(2001). 

2.2. Database of Numerical Models 

In order to gain insight into the nature of the evolu- 
tion of turbulent thermonuclear flames on large scales, 
we have constructed an extensive database of numerical 
flame models. 

Tables [0 and 2 provides a complete list of the cata- 
logued simulations. 
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TABLE 1 

Model Turbulent Thermonuclear Flames 

Parameter Definition 

N x x N y x N z model resolution a 

L lateral extent of the computational domain 

Di laminar flame speed 

Dt turbulent flame speed 

g z gravitational acceleration 

Fr Froude number, Fr = Df /g z L 

Af amplitude of the flame front perturbation in vertical direction 

A wavelength of the flame front perturbation in lateral directions 

b nominal numerical thickness of the flame front (number of computational cells) 



a Effective resolution of AMR model. 



TABLE 2 
High-Resolution non-SGS Models 3, 



Model 


N x x N y x N z 


L [cm] 


Ax 


[cm] 


Di 


; [cm s x ] 


g z [cm s~ 2 ] 


Fr x 10 4 


b 


Rl 


64 x 64 x 2560 


1.5 x 10 6 


2.34 


X 


10 4 


1 


,07 x 10 6 


1.9 x 10 9 


4.02 


4 


R2 


128 x 128 x 5120 


1.5 x 10 6 


1.17 


X 


10 4 


1, 


,07 x 10 6 


1.9 x 10 9 


4.02 


4 


R3 


256 x 256 x 10240 


1.5 x 10 6 


5.86 


X 


10 3 


1. 


,07 x 10 6 


1.9 x 10 9 


4.02 


4 


BH 


128 x 128 x 5120 


1.5 x 10 6 


1.17 


X 


10 4 


1. 


,07 x 10 6 


1.9 x 10 9 


4.02 


3 


B2 


128 x 128 x 5120 


1.5 x 10 6 


1.17 


X 


10 4 


1, 


,07 x 10 6 


1.9 x 10 9 


4.02 


6 


DH 


64 x 64 x 2560 


1.5 x 10 6 


2.34 


X 


10 4 


0. 


,54 x 10 6 


1.9 x 10 9 


1.01 


4 


D2 


64 x 64 x 2560 


1.5 x 10 6 


2.34 


X 


10 4 


2, 


,14 x 10 6 


1.9 x 10 9 


16.1 


4 


LH 


64 x 64 x 2560 


0.75 x 10 6 


2.34 


X 


10 4 


1 


,07 x 10 6 


1.9 x 10 9 


8.04 


4 


L2 


128 x 128 x 5120 


3.0 x 10 6 


1.17 


X 


10 4 


1 


,07 x 10 6 


1.9 x 10 9 


2.01 


4 


GH 


64 x 64 x 2560 


1.5 x 10 6 


2.34 


X 


10 4 


1, 


,07 x 10 6 


9.5 x 10 8 


8.04 


4 


*A f = 


: 9.96 x 10 4 cm, A = 


4.46 x 10 5 


cm. 

















Our high resolution m odels are primari ly intended to 
extend the simulations of Kho khlovl (jl995f) , verify our im- 
plementation of the flame capturing algorithm, verify the 
KSGS model, and demonstrate numerical convergence. 
Specifically, models Rl, R2 and R3 are used to investi- 
gate grid convergence. Models BH and B2 together with 
model R2 are used to demonstrate the insensitivity of the 
results to the numerical flame thicknesses, b. Models DH, 
Rl and D2 are used in a discussion of flame dynamics, in- 
cluding the self-regulation mechanism. In models LH, Rl 
and L2 we vary the size of computational domain. These 
models are analized in our study of the flame destruction 
process. Models Rl, D2, GH and L2, having various do- 
main sizes and the imposed gravitational accelerations, 
are used to verify the KSGS model. 

3. FLAME DYNAMICS 

Following jDamkohlerl ( Il939. see also Khokhlovl (|1995ft : 
iNiemever fc H illebrandt (1995)), a primary consequence 
of turbulence is to increase the flame surface area. The 
turbulent flame speed is, 



Po 



_ x dM h 
~df 



where the burned mass, M5, is normalized to the area 
of the unperturbed flame the horizontal cross-section of 
the computational domain. The turbulent flame speed is 
expected to vary in proportion to the surface area, S as 



(2) 



where Sg is th e surface area of an unperturbed flame. 
Khokhlovl (|1995l) suggested, based on examination of sim- 
ple flame surface geometries, that the evolution of the 



flame surface area can be described as a competition be- 
tween flame surface creation and destruction processes: 



where 



^-=cS-dD l S 2 . 
dt 



c = [(ei)i(ei)j + (e 2 )i(e 2 ) i ](^). 



(3) 



(4) 



Here (£1,2) i are two orthonormal vectors tangent to an 
infinitesimal surface element 5S. The first unit vector, 
ei, is orthogonal to vector normal to the flame surface, 
n = V//|V/|, and is chosen arbitrarily; the second unit 
vector, e2, is orthogonal to both ei and n and is then 
determined uniquely. 

The first term on the right-hand side of Eq. de- 
scribes surface creation by strain, while the second term 
in this equation represents the process of flame surface 
destruction due to the propagation of cusps. In steady 
state, the two terms balance one another in a statistical 
sense. Then, S oc D^ 1 follows from Eq. (|2j) implying 
that the turbulent flame speed D t is independent of the 
laminar flame speed, D\. 

Each term in Eq. (0 can be calculated using avail- 
able information from our simulations. The flame surface 
area, 5, is extracted from the data as an iso-surface of 
the progre ss variable at / = 0-5 u sing a marching-cube 
algorithm (jLorensen fc Cline1 ll987). The computational 
cells that contain flame surface (defined by / = 0.5) are 
then used to calculate the surface creation coefficient, c: 
First, a unit vector normal to the flame surface, n, is 
calculated. Then two unit vectors, ei and e p, orthog- 
onal to n and to each other are constructed (Khokhlov 
1995). It is the average value of the creation coefficient 
that is followed, J Q cdV/ J Q dV , where £1 is the volume 
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containing flame surface, i. e. all of the computational 
cells that contain flame surface. The total flame surface 
area, 5, and creation term, cS, are computed by inte- 
grating over the entire flame surface with S = J Q dS and 
cS = f Q cdS. The total destruction term is then calcu- 
lated as the difference between total creation and dS/dt, 
where dSj dt is simply calculated by finite differencing in 
time. 

3.1. Flame Surface Evolution 

We begin with an examination of the evolution of the 
model flames listed in Tables Q and 2. These results serve 
as the basis for the analysis of the actual form of the 
governing equation for flame surface evolution. First, we 
examine the evolution of models with different laminar 
flame speeds, D\. These results are closely associated 
with the self-regulation mechanism that will be discussed 
throughout this paper. Fig. [I] 



TABLE 3 

Statistics of Flame Surface in Steady State 



D, 


= 0.54x10 6 cm/s 


D, 


= 1.07 x10 6 cm/s 


D| 


= 2.14x10 6 cm/s 





D, = 0.54e6 cm/s 

- - - - D, = 1 .07e6 cm/s 
D,= 2.1 4e6 cm/s 




t, s 



7 V 



Fig. 1. — Temporal evolution of flame surface area in 2D (top) 
and 3D (bottom). Note that in 2D, the models are the counterparts 
of models DH, Rl and D2 in 3D. 



shows the evolution of the flame surface area at three 
different laminar flame speeds both in 2D and 3D. It is 



Dl xl.07 x lCT 6 cm/s (2D) S x IP- 



S' x IO- 



CS 
1. 
2. 



1.8 
1. 
0.5 



0.9 
0.5 
0.3 



Dl xl.07 x 10~ 6 cm/s (3D) S xl0~ 13 cm 



0.5. 
1. 
2. 



3.2 
2.5 
1.2 



1.2 
1 

0.5 



expected that the larger the laminar flame speed, the 
small er the flame surface. As discussed by Kho khlovl 
(1995) and shown in Sec. 13.51 the larger the laminar flame 
speed, the stronger the smoothing effect due to burning, 
thus the smaller the flame surface area. 

Even a crude inspection of Fig.^suggests that a steady 
state exists and the flame exhibits a semi-periodic behav- 
ior about this equilibrium. 4 The time averaged value and 
the standard deviation of the flame surface area in steady 
state (at t > 0.6 s) are listed in Tabled 

Both the time averaged value and the standard devia- 
tion of the flame surface area appear to scale as D^ 1 . In 
most cases considered this relation is fulfilled with accu- 
racy of ~ 10%. Following our trial time-limited high- 
resolution runs, this agreement seems improving with 
increased resolution. Both the time averaged turbulent 
flame speed and its deviation are independent of the lam- 
inar flame speed. This is in agreement with Eq. and 
is illustrated in Fig. 

Therefore, our study confirms the independence of the 
turbulent flame speed and the laminar flame speed. This 
is the essenti al property of the self-regulation mechanism 
discussed bv iKhokhlovl |l995). However, the indepen- 
dence of the deviation of turbulent flame speed versus 
laminar flame speed has not been discussed in that pre- 
vious study. 

3.2. Character of The Steady State 

As discussed in the previous section, the flame exhibits 
a semi-periodic behavior once a steady state is estab- 
lished. A natural question then is "What is the source 
of the oscillation?" It is conceivable that the oscillation 
is caused by interference between the flame and acoustic 
waves generated at the bottom wall boundary. On the 
other hand, it could be an intrinsic property of Rayleigh- 
Taylor unstable flames. In order to determine whether of 
these hypotheses is likely we perform a ID test without 
flame propagation. In this test an initial vertical down- 
ward velocity is imposed on the fluid and the propagation 
of sound waves generated at the bottom wall is exam- 
ined. We find that the sound waves quickly reach the top 
boundary and leave the computational domain. The fluid 
becomes essentially at rest once the sound wave propa- 
gates off the grid. In addition, a simulation of flame prop- 
agation in 2D with both the top and bottom boundaries 
open have been performed and the same semi-periodic 
behavior found in our fiducial models is observed. 

On the basis of these experiments, we conclude that 
flames subject to the Rayleigh-Taylor instability are in- 
trinsically oscillatory. In 2D, the rising bubbles filled 

4 The character of the steady state will be further elaborated in 
Sec. El 
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= 0.54e6 cm/s 


D, 


= 1 .07e6 cm/s 


D, 


= 2.1 4e6 cm/s 
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D 
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D 
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Fig. 2. — Evolution of the turbulent flame speed, Dt in 2D 
(left) and 3D (right) in models DH (solid line), Rl (dashed line) 
and D2 (dash-dotted line). Note that a statistically steady-state 
evolution is preceded by an extended transient. The time-averaged 
turbulent flame speed does not depend on the laminar flame speed 
once a steady-state regime is established, where the evolution is 
governed by the self-regulation of the flame surface. 

with hot burned material are analogous to a cylinder 
passing through a fluid. It is well known that this 
schematic situation typically results in th e formation o f 
a vortex street in the wake of the cylinder ( White 2003). 
Similarly, a vortex street also appears in the wake of the 
flame front in our simulations, and the periodicity of the 
structure is apparent in Fig. 03 

The production of vorticity (top panels in Fig. is 
caused by the misalignment of the density and pressure 
gradients at the sides of the flame bubbles. The pressure 
gradient tends to be oriented along the direction of grav- 
ity, whereas the density gradient is approximately normal 
to the flame surface. Vortices are formed and then shed 
from the flame surface (middle panels in Fig.|HJ). The vor- 
tices turn over, stretch and fold the flame surface near 
the tail (panels for 0.47 and 0.54 s). These stretched and 
folded surfaces then collide and annihilate. During this 
process, the heads of the bubbles wobble and their tails 
swing from side to side. Approximately one full period of 











Fig. 3. — Evolution of a turbulent flame in two-dimensions in 
statistical steady state. Shown is the distribution of vorticity gen- 
eration, Vp x Vp (top row), vorticity (middle row) and horizontal 
velocity (bottom row) at four different instants of time (from left 
to right: 0.41, 0.47, 0.54 and 0.58 s, corresponding to "A" - "D" 
marked in the first panel of Fig. EJ- Each row covers about two 
eddy turnover times ( ~ 0.2 s). Each panel is a lateral composite 
of 2.5 of our periodic domains, arranged for visualization purposes. 



this "wobbling" process is shown in Fig. [3 During each 
oscillation period, two eddies of opposite direction are 
shed from the surface of each bubble. One such turnover 
is illustrated in Fig. 0] 

The size of the eddies is comparable to the turbulent 
driving scale, L, i.e. the domain width. 

Each flame wobble is characterized by varying flame 
surface resulting in variations of the turbulent flame 
speed. For example, consider two maxima of the turbu- 
lent flame speed marked "B" and "C" in the top panel 
of Fig. [3 (2D model with D t = 1.07e6 cm/s). The flame 
structure corresponding to these two extrema is illus- 
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Fig. 4. — Eddy turnover in a turbulent flame in two-dimensions. 
Plotted is the vorticity distribution. Note the counter-clockwise 
rotating eddy (the red "blob" -like feature inside the bubble) as it 
stretches and folds the flame surface near the tail. Axes are the 
same as those in Fig. |3] 

trated in the second and the third column in Fig. 03 Note 
that, at these two instants of time, the flame surface area 
is also at its largest, as indicated by extreme length of 
the bubble's tail. 

The above picture is in qualitat ive agreement with the 
recen t two-dimensional study of IVladimirova fc Rosnerl 
(2005) obtained in the Boussinesq limit. In particular, 
the morphology of the flame and the formation of vor- 
tex streets are common featu res of the two models [see 
([Vladimirova fc Rosnerll2005l) . panels (a), (e), and (f) in 
Fig. 1]. Note that in their simulations with symmet- 
ric boundaries, no oscillatory behavior is observed. We 
consider this to be an effect of the symmetric bound- 
ary conditions. In reality, there is no imposed symmetry 
and the RT flame is physically periodic as a result of 
eddy turnover. It should be pointed out the horizontal 
reflecting boundary they sometimes term a symmetric 
boundary should not be confused with a wall bound- 
ary. Here, we observe periodic behavior also in the flame 
evolution with horizon t al wa ll boundaries. In addition, 
IVladimirova fc Rosnerl ((2Q05) surprisingly find that the 
ratio of the vertical wavelength and the domain width 
is invariant and always close to 2. A similar result is 
obtained here (Fig. EJ) for which we offer the following 
explanation: 

As mentioned above, in steady state, the largest RT 
unstable structure is always comparable to the turbu- 
lent driving scale of the system. Since the periodicity is 
caused by the turnover of these largest eddies, and there 
are two turnovers in each period, it is natural that the 
vertical wavelength is always about twice as large as the 
domain width. 

The character of the steady state has been discussed 
in the context of 2D simulations in the above. In 3D, al- 
though vortices shed from the flame surface do not merge 
into the bulk but cascade into smaller and smaller struc- 
tures, and the flame surface is not as regular (Fig. EJ, 

the overall mechanism of vorticity production near the 
flame front and the subsequent interaction between vor- 
tices and the flame surface is similar to that observed in 
2D. Regarding vorticity generation, vortex stretching is 
an important contributing mechanism of vorticty gener- 
ation uniquely present in 3D. While vortex stretching is 
a dominant source of vorticity in the bulk, its strength 



Fig. 5. — Vorticity distributions in steady state. Laminar flame 
speeds are (from left to right): 0.54e6 cm/s, 1.07e6 cm/s and 2.14e6 
cm/s corresponding to the three 2D runs shown in Fig. ^ The 
scale of the axes are kept the same in each panel. The vertical 
wavelengths are approximately the same, being about twice as large 
as the domain width. Note, again, each panel has been duplicated 
horizontally 2.5 times and axes are the same as those in Fig. 131 






Fig. 6. — Morphology of the turbulent flame surface in models 
with different laminar flame speeds at late times, (left): model 
DH, D x = 0.54 x 10 6 cm/s; (middle): model Rl, D t = 1.07 x 10 6 
cm/s; (right) model D2, D\ = 2.14 x 10 6 cm/s. Note the larger the 
laminar flame speed, the smoother and, concomitantly, the smaller 
the flame surface. The vertical spatial scale are 1 x 10 7 to 1.6 x 10 7 , 
0.7 x 10 7 to 1.35 x 10 7 and 1.2 x 10 7 to 1.8 x 10 7 cm for model 
DH, Rl and D2, respectively. 

is comparable to that of barcolinic vorticity generation 
near the flame surface. The latter mechanism is common 
to 2D and 3D and is the original source of vorticity in 
our problem. 

In addition, the same semi-periodic behavior is also 
observed in 3D flame evolution, as can be seen in Fig. [3 
Just as in 2D, this semi-periodic behavior is caused by 
the turnover of largest eddies. Despite the complicated 
structure of a 3D flame, the periodically layered structure 
can still be discerned in a 2D slice of the 3D flow by 
coarsening the grid by a factor of 16. Such slices are 
shown in Fig. El 

To actually extract the period of the temporal evolu- 
tion of the flame, the run shown in (Fig. |2J) was ana- 
lyzed for t > 0.6 s in both 2D and 3D. Figure |S| shows 
the results of the anal ysis obtained for 3D with the help 
of the CLEAN method Roberts, Lehard, DrehsJ|1987|) 
for power spectrum analysis, together with the results 
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Fig. 7. Turbulent flame at ~0.83 s for model Rl. The 3D 
flame surface and the position of the 2D slice are shown in the 
first panel. Only that part of the domain that contains the flame 
(vertical spatial scale from 1 x 10 7 to 1.6 x 10 7 cm) is shown in 
this first panel. The second panel plots the horizontal velocity 
field (component in the plane) in the 2D slice. The 3rd and 4th 
panel show the horizontal velocity field and the vorticity (normal 
to the plane) on the coarsened grid. Note the trace of the layered 
structure in these two panels. (N.B. In contrast to the 3D plot, the 
slices cover the entire domain height from to 3 x 10 7 cm.) 

obtained by computing the auto-correlation of the time 
series D t (at t > 0.6 s) in Fig. , i.e. the correlation of 
D t (t) against a time-shifted (by r) version of it. Fig. [HI 

shows that both methods indicate the presence of ex- 
cess power around 0.25. In particular, there is a clear 
accumulation of power for frequencies < 10 s _1 in the 
power spectrum calculation. The first and second peaks 
of the auto-correlation function are separated by ~ 0.2 s, 
in excellent agreement with presence of power at a fre- 
quency ~ 5 8 _1 in the clean result. Fig. El 

shows the corresponding clean results for three se- 
lected 2D models. Since, as discussed in Sec. 13.11 there 
are two eddy turnovers during one period, we expect 
the power to be concentrated at low frequencies corre- 
sponding to one (half period) and two (full period) eddy 
turnover times. This expectation is confirmed, as can be 
seen in Fig. [HI 

This period determination is consistent with the time 
for the flame "brush" to travel the distance between the 
vortex streets left in the bubble wake. In the particu- 
lar case shown in Fig. 03 the distance between two con- 
secutive, countervailing vortex streets is ~ 2.7 x 10 6 cm 
(this can also be measured in the bottom panels of Fig. 01 
more easily), while the flame sweeps through the compu- 
tational domain with total speed ^ 1.5 x 10 7 cm/s. 5 

Our results in this section lead to the conclusion that 
the semi-periodic behavior we observe is caused by the 
turnover of the largest eddies present in the simulation. 
This behavior is clearly visible in 2D simulations. In 
the supernova setting, a steady state realized in our 
model flame study is, strictly speaking, not present due 
to supernova expansion that makes both gravity and 
pre-flame conditions time-dependent. Therefore, super- 

5 The total speed is the sum of the turbulent flame speed, ~ 
0.8 x 10 6 cm/s (Fig.0, and the mean fluid velocity, ~ 0.7 x 10 6 
cm/s. 




T, S 



Fig. 8. — Time scales for the evolution of turbulent flames. 
The CLEAN power spectrum (left) and the auto-correlation function 
R(t) (right) are plotted for the 3D time series data shown in Fig.|2| 
In the first panel, the power peaks at a frequency of « 5 s -1 . This 
is consistent with the separation of the first and the second peak 
of the auto-correlation by « 0.2 s. 

nova flame never achieves steady state but rather con- 
stantly struggles to approach that (now) time-dependent 
condition. One may expect that in that case a semi- 
periodic behavior exhibited by model flames might be 
significantly weaker or perhaps even completely sup- 
pressed (depending on whether flame evolution proceeds 
on timescales short compared to that of the background 
expansion). 

3.3. Governing Equation for Flame Surface Area 
Evolution 

As previously mentioned, Khok Movl (|1995l) suggested 
that the flame surface evolution is governed by Eq. ©. 
The first term on the right-hand side of this equation 
describes flame surface creation via strain, while the sec- 
ond term describes fl ame surfac e destr uction due to the 
propagation of cusps. Khokhlov (1995) motivated a gen- 
eral S 2 dependence of the flame surface destruction term 
by considering elemental geometries such as intersecting 
spheres and parallel planes. 

To determine whether the flame surface destruction 
term in (0) is indeed proportional to S 2 let us consider 
the following generalized equation: 

^=aS-bS n , (5) 

where a = c and b = dD\ are arbitary constants for 
the moment, and n > 0. A steady state (equilibrium) 



8 



Zhang et al 





Fig. 9. — Evolution of the turbulent flame speed, Dt in 2D 
in steady state (left panel) and the corresponding CLEAN power 
spectrum (right panel). Note that offsets of 0.9 x 10 7 and 1.8 x 
10 7 cm/s have been added to the slowest and the second slowest 
laminar flame speed cases, respectively. 

solution, S e , to (|5J) can be found by equating the creation 
(first) and destruction (second) term to produce, 



(6) 



Now, let us add a constant perturbation Aa to the sur- 
face creation coefficient, a, to obtain a new equilibrium 
solution S' as, 



Aa 



, (7) 



Therefore, the difference between the two equilibrium so- 
lutions is 



AS =S' e -S e 

= (f)^[(l + ^)^-l]- 

In the limit of < 1, Eq. © yields 

Aa / a \ tt~=t 



(8) 



AS 



(n-1) 



a (5) 



(9) 



/From this one can see that for a given a, b and n > 1, 
the magnitude of AS is proportional to Aa. Also, for a 

given a, Aa and n > 1, AS is proportional to (|) ™ _1 • 

Let us now consider a time-dependent perturbation, 
Aa(t), and ask if the above properties still obtain under 
this generalization. Assume 



Aa(t) = Aao sin(a;t) . 



(10) 



where Aao and co are the amplitude and frequency of the 
perturbation, respectively. We have obtained numerical 
solutions for a = 2 and Aao = 0.2. Three different val- 
ues of n and b (equivalently, the laminar flame speed) 
were examined to determine the dependence of the equi- 
librium solution, S e , and the deviation from it, AS, on n 
and b. Fig.lTUlshows these solutions. The dependences of 
S e and AS on n and b described by Eqs. © and (0 re- 
main the same under the generalization of perturbation. 
Specifically, for a given a and Aao, we always obtain 



(!)* 



AS- 



Aao f a \ n"=T 



(n 



Da \b) 



(11) 



(12) 



With these two functional dependences firmly estab- 
lished, we are left with determining the value of n in 
© to completely specify the governing equation of the 
flame surface area evolution. 

As was demonstrated in Sec. 13.11 both the time- 
averaged value and the standard deviation of the flame 
surface area is proportional to 1/Di. Eqs. (|TT1 and (|T^1 
then imply that n, in Eq. © is indeed equal to 2. (A new 
subgrid model can then be developed based on Eq. Q 
of the flame surface area evolution. Details of the new 
model will be discussed in a forthcoming paper.) There- 
fore, we restrict our immediate discussion to the case 
n = 2. Note that we assume the creation coefficient, a, 
and its deviation, Aao, are constant for different laminar 
flame speed, D^s. We will provide evidence that this is 
a reasonable assumption in the following section. 

In case n = 2, Eqs. (H]) and O yield 



AS 



a 
~~b' 

Aa 



(13) 



(14) 



The amplitude of the variation of the flame surface 
area, S, in the steady state is proportional to the am- 
plitude of the variation of the perturbation, Aao- In- 
tuitively, there should exist a relaxation time scale, r r , 
associated with the response of S to a perturbation in 
the creation coefficient, a. For a = const, the solution of 
Eq. © is 



SoS e 



S e e~ at + S (l 
and can be also written as 

a-bS S e -S S, 



So 



bS 



S 



So 



(15) 



(16) 



where So = S(t = 0). Therefore, the relaxation time 
scale r r is determined solely by the value of the coefficient 
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: 2, b = 1 
= 2, b = 2 
= 2, b = 4 



300 
t, S 




t, S 



a=2, b = 1 
a=2, b = 2 
a=2, b = 4 



t, S 



Fig. 10. — Solution of Eq. J3 with different b and n's. Top: 
n = 1.5; middle: n = 2; Bottom: n = 3. Creation = a(l. + 
0.1sin(o;t))S for all cases. 

a itself. For specificity, let us assume r r = 1/a. Note 
that this relaxation time scale does not depend on b (or, 
equivalently, the laminar flame speed). We can see from 
Fig. El that, for the perturbation defined in Eq. (|10[). 



there is no phase difference between cases with different 
laminar flame speed. The difference in laminar flame 
speed affects only the steady state value of S and the 
deviation from it. 



0.1, T = 62.8 s 

0. 5, T = 12.6 s 

1, T = 6.28 s 

2, T = 3.14 s 
10, T = 0.628 s 
100, T = 0.063 s 




Fig. 11. — Solution of Eq. J5j for several values of uo. a =2, b =1, 
and Aao = 0.2 for all cases. Note the percentage variation in flame 
surface area, S, is close to 10% only when the period of excitation, 
T, is small compared to the relaxation time scale r r = 1/a = 0.5 
s. 



Given that it takes some time for the flame surface 
to respond to a change in the flow conditions, we might 
expect the variation of the flame surface area to be differ- 
ent when the frequency of the excitation, cj, in Eq. 
changes. To examine the response of the flame surface 
area to excitations of different frequencies, a numeri- 
cal experiment was performed with a = 2, b = 1 and 
Aao = 0.1a = 0.2. The frequency of the excitation is var- 
ied from 0.1 to 100 s" 1 . Combining Eqs. lO and 
we obtain 

AS _ Aa 

S e a 

i.e., the percent variation in the flame surface area is 
equal to the percent variation of the creation coefficient. 
It is clear from Fig. El that only when the frequency of 
excitation is low (uu < 1) is the amplitude of the variation 
in S as large as the percentage variation of the creation 
coefficient, i.e. ~10%. When the frequency of excitation 
is large (u > 1), the flame surface becomes less respon- 
sive and the percentage by which the flame surface area 
varies decreases to below 10%. The relaxation time scale 
in this case is r r = 1/a = 0.5 s. The period of excita- 
tion, T, is 6.28 s for an excitation frequency of 1 s _1 . 
Therefore, the period of excitation needs to be about ten 
times larger than the relaxation time scale for the flame 
surface to "fully" respond to the excitation. 

3.4. Flame Surface Creation and Destruction 

After determining the form of the governing equation 
of the flame surface area evolution, let us now attempt to 
determine the parameters in the equation for the actual 
physical system. 
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As introduced at the beginning of Sec. 03 the evolution 
of the creation coefficient, c, is followed in terms of its av- 
erage value, J n cdV/ f n dV where Q is the total volume 
containing flame surface, i. e. all of the computational 
cells that contain flame surface. Fig. fT21 



0.54e6 cm/s 




20 - 



D, 


= 0.54 x 1 6 cm/s 


D, 


= 1.07 x10 6 cm/s 


D, 


= 2.14 x10 6 cm/s 



o ' 

0.2 0.4 0.6 0.8 1 

t, S 



Fig. 12. — The temporal evolution of the surface creation in 2D 
(left) and 3D (right) in models DH, Rl and D2 and their counter- 
parts in 2D. Note the similar time averaged value of c at different 
laminar flame speeds. 

shows the evolution of the average creation coefficient 
for models DH, Rl and D2. One can see that c varies 
around 30 s _1 in 2D and 60 s _1 in 3D at different lami- 
nar flame speeds. 6 As discussed in the previous section, 
the reciprocal of the creation coefficient, a = c, is the re- 
laxation time scale (r r ) of flame response to a change in 
the flow conditions. The relaxation time scale is « 0.033 
and 0.017 s, in these models for 2D and 3D, respectively, 
while the RT time scale, trt = y/L/Ag, is 0.05 s. In 
addition, c « 50 s _1 in a model with nominal gravity 
while about twice smaller than that in the model with 
gravity 4 times lower (see Fig. Et - 

6 Note the lower creation coefficient in the D2 run, where the 
flow field is not quite as turbulent as the other runs (cf Figs. |5| 
and|6|. 




g = g /4 



o ' 

0.1 0.2 0.3 0.4 0.5 

t, S 



Fig. 13. — The temporal evolution of the surface creation coeffi- 
cient in 3D models with nominal (solid) and 4 times lower gravity 
(dashed). 

Therefore, the reciprocal of the creation coefficient 
seems to also be a measure of the RT time scale. This 
is actually the physical explanation for the independence 
of the relaxation time scale from the laminar flame speed 
discussed in the previous section. 

It is also interesting to note that the creation coeffi- 
cient in 3D is about twice as large as that in 2D. It can 
be shown analytically that the creation coefficient for a 
uniformly expanding sphere is twice as large as that of an 
expanding cylinder. It appears this effect of dimension- 
ality can be generalized to more generic shapes. The ad- 
ditional potential for surface area creation afforded by an 
additional dimension also serves to explain why our de- 
termined turbulent flame speeds are roughly twice larger 
in 3D than in 2D (c/Figs. EJ). 

The numerical value of the surface destruction coeffi- 
cient, can be estimated from the simulation data using 
Eq. (O- We obtain the destruction term by subtracting 
the rate of change of the flame surface area from the 
flame surface creation term, d is then obtained by divid- 
ing the destruction term by D\S 2 . Dimensional analysis 
indicates that the coefficient d must have the dimensions 
of specific volume. It is therefore reasonable to expect 
that d scales as 1/L 3 , where L is the lateral extent of the 
computational domain. To verify this result, we calcu- 
late three models of turbulent flame evolution in boxes 
of different lateral dimension, each model having lateral 
extent twice that of its predecessor (models LH, Rl and 
L2). The evolution of the destruction coefficient in those 
runs is shown in Fig. EI 

which clearly indicates that the destruction coefficient 
scales with the inverse of the simulation volume, i.e. d oc 
1/L 3 . Given that the Ray leigh- Taylor bubble volume is 
comparable to L 3 at late times (Fig. H5|) . one can state 
that the destruction coefficient scales with the inverse of 
the characteristic Ray leigh- Taylor bubble volume. 

The following interpretation of destruction coefficient 
scaling with the inverse of the Ray leigh- Taylor bubble 
volume can be offered: 

Consider a box of size L per dimension filled with seg- 
ments of flame surface. Let us assume that the segments 
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Fig. 14. — Comparison of the time evolution of dL 3 in models 
LH, Rl and L2, where d is the destruction coefficient and L is 
the lateral extent of the computational domain. Note that the 
destruction coefficient scales roughly as 1/L 3 . 





Fig. 15. — Morphology of the turbulent flame surface in models 
with different domain sizes at late times, (left) model Rl: L = 
1.5 x 10 6 cm; (right) model LH: L = 0.75 x 10 6 cm. Note that 
at late times the bubble size, becomes comparable to the size of 
computational domain, L. 

are parallel to each other and are separated by distance 
A. Moreover, let us assume that the neighboring flame 
segments are moving towards each other with a constant 
speed, D\. Then the surfaces are destroyed in collisions 
at the rate 



S 2 



A/A SA/Di ' 
In this case, the destruction coefficient is 



d = 



1 

5A 



which can been seen as a measure of how much flame 
surface can be packed into a given volume. The term 
SA can be identified with the bubble volume. 



The fact that the destruction coefficient is inversely 
proportional to the bubble volume, combined with di- 
mensional analysis of the surface area change rate equa- 
tion strengthens our conclusion that the surface destruc- 
tion term is indeed prop ortional to S 2 , in agreement with 
the original proposal of Khok hlovl (|l995l) . 

3.5. Froude Number and the Self- Regulation Mechanism 

The flame surface creation and destruction processes 
reflect the interplay between the Ray leigh- Taylor insta- 
bility and the tendency of the burning process to smooth 
the flame's surface. This interplay has been clearly il- 
lustrated in Fig. ED where we show the flame surface 
in three models obtained with different laminar flame 
speeds. Model DH (left panel in Fig. |5J) has the lowest 
laminar flame speed and develops a flame surface which is 
substantially more convoluted than in the remaining two 
cases where the laminar flame speed is higher (models Rl 
and D2, middle and right panel in Fig. respectively). 
The flame surface is relatively smooth in high laminar 
speed model D2. This behavior is well characterized by 
the Froude number, 

= _L 

gL 

which is the only parameter required to completely char- 
acterize the behavior of our model flames provided the 
density is constant. In models with higher Froude 
number the flame surface smoothing due to burning is 
stronger, resulting in a less convoluted flame surface. 
Conversely, in models with lower Froude number, the 
flame surface is strongly convoluted. The burning rate, 
however, remains controlled by the self-regulation mech- 
anism and is independent of Froude number (i.e. inde- 
pendent of the laminar flame speed). 

^From the point of view of flame surface creation pro- 
cesses, systems characterized by smaller Froude number 
are expected to have similar surface creation coefficients, 
independent of D\ as long as the RT time scale are simi- 
lar (as seen in Fig.^J). However, the creation coefficient 
will become sensitive to the laminar flame speed once the 
Froude number is such that the flame speedup factor, i.e. 
the ratio of turbulent flame speed to the laminar flame 
speed, becomes less than ~10 for 3D. 

3.6. Verification of the KSGS model 
3.6.1. Numerical Convergence 

To examine the convergence behavior of our numerical 
flame models, we perform simulations at three different 
spatial resolutions (models Rl, R2 and R3). Fig. fTTH 

shows the evolution of the turbulent flame speed, D tl 
in these three models. At sufficiently late times the flame 
speed exhibits variations of ~ 15% around a mean value 
of & 1.6 x 10 T cm/s, and the detailed behavior is remark- 
ably similar in all three cases. (Note that, due to limited 
computing resources, the highest resolution calculation 
was stopped at t w 0.5 s. Nevertheless, even in this case 
the flame speed increase saturates at a similar level as in 
the two lower resolution cases.) We also find that at any 
given instant of time the total burned mass is similar in 
all three models. 

Due to the inviscid nature of the Euler equations, the 
flow structure will depend on numerical resolution to a 
degree. This is illustrated in Fig. IT7I 
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6 (model B2). Figure d 




Fig. 16. — Convergence study for the turbulent flame models. 
Shown is the temporal evolution of the turbulent flame speed, Dt : 
in models Rl, R2 and R3. 




Fig. 17. — Dependence of the morphology of the turbulent flame 
surface on numerical resolution. The surface of the flame is shown 
at t 0.2s in model Rl (left panel, Ax = 2.34 x 10 4 cm), R2 
(middle panel, Ax = 1.17 x 10 4 cm) and R3 (right panel, Ax = 
5.86 x10 s cm). 

where we show the morphology of the flame surface in 
models Rl, R2, and R3 obtained at increasing resolution 
(but having fixed Froude number). One can notice a dra- 
matic increase in the amount of small-scale structure as 
the resolution is increased. However, once the large scale 
characteristics of the flame are captured by the numerical 
model, the small scale structure emerging with increased 
numerical resolution does not change the global proper- 
ties of the evolution, as seen in the convergence of the 
turbulent flame speed shown in Fig. El 

Our results might also be expected to depend on cer- 
tain technical aspects of the simulations. For example, 
due to the finite thickness of the numerical flame front it 
is desirable to study the dependence of our results on this 
numerical parameter. For this purpose we obtained mod- 
els BH and B2, in which we adjust the nominal numerical 
flame thickness to, respectively, 0.75 and 1.5 times the 
standard value used in model Rl. This is accomplished 
by changing the parameter b from 4 to 3 (model BH) or 




Fig. 18. — Dependence of the turbulent flame speed on the 
nominal numerical thickness of the flame front, 6, in models BH, R2 
and B2. Note that in all cases the turbulent flame speed oscillates 
around a similar average value. 

shows that at late times the turbulent flame speed in 
all three models oscillates around a similar mean value 
with an amplitude entirely consistent with that observed 
in the resolution study discussed above. 

Since our fiducial simulations make use of adaptive dis- 
cretization, it is necessary to verify that local nonunifor- 
mity of the computational grid does not markedly influ- 
ence the evolution of the system. Figure IT9l 




Fig. 19. — Verification of the adaptive mesh refinement model 
of the turbulent flame. Temporal evolution of the turbulent flame 
speed in the AMR model is plotted along that obtained in the 
simulation with uniform grid. The two models produce statistically 
similar results. 

shows the evolution of the turbulent flame speed in 
model Rl and its counterpart calculation on a uniform 
mesh. The two models perfectly agree at very early times 
(t < 0.1s) when the flame evolution remains confined 
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to the region initially refine. The two models begin to 
differ during the initial transient phase when the flame 
structure becomes progressively turbulent and the mesh 
in model Rl begins to follow the evolving flame structure. 
However, in steady state, the flame speeds in the two 
models saturate at similar levels. Also, we find that the 
power spectrum of the time variation of the flame surface 
area obtained in the model calculated on a uniform grid 
shows the same salient features (location of peak power, 
concentration at long timescales) as the AMR-enabled 
model. 

3.6.2. Examining the Self- Regulation Mechanism 

One of the main goals of this stud y is to verify th e 
KSGS subgrid-scale model proposed bv lKhokhlovl Jl995). 
The KSGS model considers turbulence driven by the 
Ray leigh- Taylor instability on scales comparable to the 
grid resolution. In this model, the instantaneous turbu- 
lent flame speed, D tl is 

A = 0.5v^Z, (17) 

where A = (po — pi)/ (po + pi) is the Atwood number, po 
and pi are the densities ahead and behind the flame front, 
respectively, g is the gravitational acceleration, and L is 
the turbulent driving scale, usually set to one or two com- 
putational cell widths i n large-scale simula tions making 
use of the KSGS model ([Gamezo et al.ll2003l) . A number 
of models were calculated to verify the applicability of 
Eq. (|T7|) for different g and/or L . The results are shown 
in Fig. EH 




Fig. 20. — Verification of the KSGS model of RT-driven turbulent 
flames. The time evolution of the turbulent flame speed is shown 
for models Rl (reference model), D2 (twice laminar flame speed), 
GH (half gravitational acceleration), and L2 (twice larger domain 
size). Time averaged turbulent flame speeds in the stea dy-state 
regime are within 10% of the values predicted by Eq. 1171 . 

The time-averaged steady-state turbulent flame speeds 
are entirely consistent with the predictions of Eq. (|T7)) . 

Another aspect of the verification of the KSGS model is 
to confirm the independence of the self-regulation mech- 
anism from details of the flame microphysics, i.e. on the 
laminar flame speed. To verify this, models DH, Rl, and 



D2 were obtained for the same initial conditions and res- 
olutions, but using laminar flame speeds differing by a 
factor of 2. The long-term evolutions of the turbulent 
flame speeds in these models are shown in Fig. |2j which 
clearly illustrates the expected independence of the tur- 
bulent and laminar flames speeds. 

4. SUMMARY 

We present the results of an extensive study of the 
evolution of thermonuclear flames in periodic domains 
under the influence of a constant gravitational field. The 
flame surface is Ray leigh- Taylor unstable, leading to a 
growth of initial perturbations. We follow the flame's 
development through an extended initial transient phase 
and find a well-defined steady-state regime. We verify 
that the evolution of RT-unstable flames is self-regulated, 
i. e. the turbulent flame surface evolution in steady state 
is independent of the details of the incorporated micro- 
physics. 

The evolution of the flame surface was examined. The 
functional dependence of flame surface destruction sug- 
gested by (Khokhlov 1995) is proven and the properties 
of the governing equation of the surface evolution are 
discussed. The flame surface creation coefficient is found 
to be closely related to the RT time scale and the flame 
surface destruction strength is found to vary as l/I/ 3 , 
where L is the turbulent driving scale. The flame surface 
creation and destruction processes reflect the interplay 
between the Rayleigh-Taylor instability and the flame's 
tendency to smooth its surface. We find this relationship 
can be well characterized by the Froude number. 

We find that the turbulent flame speed obeys a \J AgL 
scaling law for Froude numbers <C 1. Also, we are able 
to identify semi-periodic temporal variations in the flame 
evolution in steady state and associate them with the 
turnover time of the largest eddies; no significant amount 
of power on shorter time scales is found, indicating the 
evolution is governed by large-scale flow characteristics. 

We discussed in some detail the mechanism of vorticity 
generation, shedding, and transport in the vicinity of the 
flame surface. We find that the interaction between in- 
duced vorticity and the flame surface plays an important 
role in the flame evolution, leading to bubble "wobbling" 
on timescales comparable to the eddy turnover time. As 
the bubbles rise, vorticity is advected downstream into 
the bulk and forms an homogenized, mixed turbulent 
layer. It is in this layer where most of the flame surface 
destruction takes place. Although there exisst notice- 
able differences between flame dynamics in 2D and 3D, 
the vorticity production and advection, as well as the in- 
teraction between vortices and the flame surface appear 
quite generic. 

There are several possible future directions of research 
emerging from our study. Our convergence study indi- 
cates that the self-regulation mechanism is relatively in- 
sensitive to the numerical resolution, provided the driv- 
ing scales are adequately resolved. However, at the high- 
est tested resolutions, we exhaust our current computing 
resources and are unable to investigate evolution with the 
Gibson scale fully resolved. Such investigation, bridging 
large and small scales in one model, is highly desirable. 
A full 3D implementation and tests including large-scale 
Type la simulations with a new subgrid model will also 
be conducted. 
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APPENDIX 



In order to verify the correctness and applicability of the CLEAN algorithm implementation, we have generated 
and analyzed a set of trial time series. The first test time series consists of single frequency signal with v = 19.5 s _1 . 
Figure Q 
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Fig. 1. — CLEAN power spectrum of trial time series with single frequency without (solid) and with (dotted) frequecy drift. See text for 
details. 



shows that this single frequency is clearly identified by CLEAN (solid line). We can extend this trivial example 
to consider the effect of frequency drift on our ability to identify the spectrum. In a second test the time variable 
is transformed using a simple power-law, £(°- 9 +°- 2x (^Ae)) ? where t e is the final time. The resulting power spectrum is 
shown in Fig. (dotted line). The power now occupies a much wider range of frequencies, consistent with the assumed 
drift of the base frequency. We posit that such a drift around a base frequency occurs in the complex interaction 
between the RT-unstable flame with its own induced turbulent velocity field. This leads to a substantial leakage of 
power to frequencies neighboring the characteristic RT (turnover) frequency. 

In another verification test, the trial signal consists of multiple discrete frequencies superposed on a white noise 
signal: 

N 

s(t) = A r (t) £ Ai[l + n(t)] (508(71-14(2* + l(T 2 r^))), 
i = 1 

where N = 4, r| ,i} and r' are time-dependent random numbers between -1 to 1, (Ai,Vi) (chosen to be are (0.1,4), 
(0.1,5), (0.3,19.5), (0.1,38)) are the amplitudes and matched frequencies of the discrete signals, respectively, and the 
amplitude of the white noise is Ao = 0.1. The time series consisted of 100 equally spaced points in time between 
and 1. Despite the fact that (a) some frequencies are closely spaced and/or aliased, (b) the amplitude of white noise is 
comparable to that of the composite signals, and (c) small random perturbations have been added in the time domain, 
most all the informational content in the signal is identified by the CLEAN method (Fig. EJ). 

Only in the case of the two lowest closely-spaced frequencies (y =4,5) does the method fail to clearly delineate the 
separate signals. Even in this case the power peak is relatively broad, suggesting the presence of distinct signals in the 
region. 
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Fig. 2. — CLEAN power spectrum for signal with multiple, closely-spaced, and aliased frequencies in the presence of significant white 
noise. Note that all major component frequencies are clearly identified with the exception of the two low frequency signals (power peak 
around v = 5 s _1 ). 
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